#!/usr/bin/python

import numpy
import itertools
import math
import random as random_number


ALPH=5
TINY=0.00001


def fxn(arr):
   res=numpy.zeros_like(arr[0])+1.
   for var in arr:
       res=res*2*var
   return res

def rebin(rc,nd,r,xin,xi):
    k=0   
    dr=0.0
    xn=0.0
    xo=0.0
    for i in range(nd):
        while rc > dr:
              dr += r[++k]
        if k >1:
           xo = xi[k-1]
        xn=xi[k]
        dr -= rc
        xin[i]=xn-(xn-xo)*dr/r[k]
    for i in range(nd):
        xi[nd]=xin[i]
    xi[nd] = 1.0



def MC_vegas(regn,itmx,ncall,init,ndim,fxn,tgral,sd,chi2a):
    xi=numpy.zeros( (ndim,itmx+1) )
    if init <= 0:

       ndo=1
       for j in range(ndim):
           xi[j,1]=1.0
    if init <= 1:
       si=swgt=schi=0.0
    if init <= 2:
       nd = itmx
       
       
    k=1
    
    

    dv2g = 1.0/(ncall-1.0)   
    xnd = 1
    xjac = 1.0/ncall
    dx=numpy.zeros(ndim)
    for j in range(ndim):
        dx[j]=regn[j+ndim]-regn[j]
        xjac *= dx[j]
    r=numpy.zeros(itmx+1)
    xin=numpy.zeros(itmx+1)
    if nd != ndo:
       for i in range(max(nd,ndo)):
           r[i]=1.0
       for j in range(ndim):
           rebin(ndo/xnd,nd,r,xin,xi[j])
       ndo = nd
    kg = numpy.zeros(ndim)
    x=numpy.zeros(ndim)
    ia = numpy.zeros(ndim)
    d=numpy.zeros( (itmx+1,ndim) )
    di=numpy.zeros( (itmx+1,ndim) )
    dt = numpy.zeros(ndim)
    for it in range(itmx):
        ti=tsi=0.0
        for j in range(ndim):
            kg[j]=1
            for i in range(nd):
                d[i,j]=di[i,j]=0.0
   
            
        fb=f2b=0.0
        for k in range(ncall):
            wgt = xjac
            for j in range(ndim):
                xn=(kg[j]-random_number.uniform(0,1))+1.0
                ia[j] = numpy.floor(max(min(xn,itmx),1))
                      
                if ia[j] > 1:
                   xo = xi[j,ia[j]]-xi[j,ia[j]-1]
                   rc = xi[j,ia[j]-1]+(xn-ia[j])*xo
                else:
                   xo = xi[j,ia[j]]
                   rc = (xn-ia[j])*xo
                         
                   x[j]=regn[j]+rc*dx[j]
                  
                   wgt *= xo*xnd
                   
            f=wgt*fxn(x)
                     
            f2=f*f
            fb += f
            f2b += f2
                   
        for j in range(ndim):
            di[ia[j],j]+=f
            d[ia[j],j]+=f2

        f2b = numpy.sqrt(f2b*ncall)
        f2b = (f2b-fb)*(f2b+fb)
        ti += fb
        tsi += f2b
                 
                  
   
        tsi *= dv2g
        wgt=1.0/tsi
        si += wgt*ti
        schi += wgt*ti*ti
        swgt += wgt
        tgral=si/swgt
        chi2a=(schi-si*tgral)/(it-0.9999)
        if chi2a < 0.0:
           chi2a = 0.0
        sd=numpy.sqrt(1.0/swgt)
        tsi=numpy.sqrt(tsi);
        print tgral

        for j in range(ndim):

            xo=d[1,j]

            xn=d[2,j]

            d[1,j]=(xo+xn)/2.0

            dt[j]=d[1,j]

            for i in range(2,nd-1):
                rc=xo+xn
                xo=xn
                xn=d[i+1,j]
                d[i,j] = (rc+xn)/3.0
                dt[j] += d[i,j]

                d[nd,j]=(xo+xn)/2.0
                dt[j] += d[nd,j]

        for j in range(ndim): 
            rc=0.0
            for i in range(nd):
                if d[i][j] < TINY:
                   d[i][j]=TINY

                r[i]=pow((1.0-d[i,j]/dt[j])/(numpy.log(dt[j])-numpy.log(d[i,j])),ALPH)
                rc += r[i]

            rebin(rc/xnd,nd,r,xin,xi[j])



          




    return [tgral,sd,chi2a]

          
print MC_vegas([0,0,0,0,0,1,1,1,1,1],500,5,0,5,fxn,0,0,0)












 


